home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / scale_geom / filt.c next >
Encoding:
C/C++ Source or Header  |  1991-09-17  |  11.2 KB  |  477 lines
  1. /*
  2.  * filt: package of 1-d signal filters, both FIR and IIR
  3.  *
  4.  * Paul Heckbert, ph@miro.berkeley.edu    23 Oct 1986, 10 Sept 1988
  5.  *
  6.  * Copyright (c) 1989  Paul S. Heckbert
  7.  * This source may be used for peaceful, nonprofit purposes only, unless
  8.  * under licence from the author. This notice should remain in the source.
  9.  *
  10.  *
  11.  * Documentation:
  12.  *  To use these routines,
  13.  *    #include <filt.h>
  14.  *  Then call filt_find to select the desired filter, e.g.
  15.  *    Filt *f;
  16.  *    f = filt_find("catrom");
  17.  *  This filter function (impulse response) is then evaluated by calling
  18.  *  filt_func(f, x).  Each filter is nonzero between -f->supp and f->supp.
  19.  *  Typically one will use the filter something like this:
  20.  *    double phase, x, weight;
  21.  *    for (x=phase-f->supp; x<f->supp; x+=deltax) {
  22.  *        weight = filt_func(f, x);
  23.  *        # do something with weight
  24.  *    }
  25.  *
  26.  *  Example of windowing an IIR filter:
  27.  *    f = filt_find("sinc");
  28.  *    # if you wanted to vary sinc width, set f->supp = desired_width here
  29.  *    # e.g. f->supp = 6;
  30.  *    f = filt_window(f, "blackman");
  31.  *  You can then use f as above.
  32.  */
  33.  
  34. /*
  35.  * references:
  36.  *
  37.  * A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975
  38.  *
  39.  * R.W. Hamming, Digital Filters, Prentice-Hall, Englewood Cliffs, NJ, 1983
  40.  *
  41.  * W.K. Pratt, Digital Image Processing, John Wiley and Sons, 1978
  42.  *
  43.  * H.S. Hou, H.C. Andrews, "Cubic Splines for Image Interpolation and
  44.  *    Digital Filtering", IEEE Trans. Acoustics, Speech, and Signal Proc.,
  45.  *    vol. ASSP-26, no. 6, Dec. 1978, pp. 508-517
  46.  */
  47.  
  48. static char rcsid[] = "$Header: filt.c,v 2.2 88/12/30 21:32:17 ph Locked $";
  49. #include <math.h>
  50.  
  51. #include "simple.h"
  52. #include "filt.h"
  53.  
  54. #define EPSILON 1e-7
  55.  
  56. typedef struct {    /* data for parameterized Mitchell filter */
  57.     double p0, p2, p3;
  58.     double q0, q1, q2, q3;
  59. } mitchell_data;
  60.  
  61. typedef struct {    /* data for parameterized Kaiser window */
  62.     double a;        /* = w*(N-1)/2 in Oppenheim&Schafer notation */
  63.     double i0a;
  64.     /*
  65.      * typically 4<a<9
  66.      * param a trades off main lobe width (sharpness)
  67.      * for side lobe amplitude (ringing)
  68.      */
  69. } kaiser_data;
  70.  
  71. typedef struct {    /* data for windowed function compound filter */
  72.     Filt *filter;
  73.     Filt *window;
  74. } window_data;
  75.  
  76. void mitchell_init(), mitchell_print();
  77. void kaiser_init(), kaiser_print();
  78. void window_print();
  79. double window_func();
  80. static mitchell_data md;
  81. static kaiser_data kd;
  82.  
  83. #define NFILTMAX 30
  84. static int nfilt = 0;
  85.  
  86. /*
  87.  * note: for the IIR (infinite impulse response) filters,
  88.  * gaussian, sinc, etc, the support values given are arbitrary,
  89.  * convenient cutoff points, while for the FIR (finite impulse response)
  90.  * filters the support is finite and absolute.
  91.  */
  92.  
  93. static Filt filt[NFILTMAX] = {
  94. /*  NAME    FUNC        SUPP    BLUR WIN CARD UNIT  OPT... */
  95.    {"point",    filt_box,    0.0,    1.0,  0,  1,  1},
  96.    {"box",    filt_box,       0.5,    1.0,  0,  1,  1},
  97.    {"triangle",    filt_triangle,  1.0,    1.0,  0,  1,  1},
  98.    {"quadratic",filt_quadratic, 1.5,    1.0,  0,  0,  1},
  99.    {"cubic",    filt_cubic,     2.0,    1.0,  0,  0,  1},
  100.  
  101.    {"catrom",    filt_catrom,    2.0,    1.0,  0,  1,  0},
  102.    {"mitchell",    filt_mitchell,    2.0,    1.0,  0,  0,  0,
  103.     mitchell_init, mitchell_print, (char *)&md},
  104.  
  105.    {"gaussian",    filt_gaussian,  1.25,    1.0,  0,  0,  1},
  106.    {"sinc",    filt_sinc,      4.0,    1.0,  1,  1,  0},
  107.    {"bessel",    filt_bessel,    3.2383,    1.0,  1,  0,  0},
  108.  
  109.    {"hanning",    filt_hanning,    1.0,    1.0,  0,  1,  1},
  110.    {"hamming",    filt_hamming,    1.0,    1.0,  0,  1,  1},
  111.    {"blackman",    filt_blackman,    1.0,    1.0,  0,  1,  1},
  112.    {"kaiser",    filt_kaiser,    1.0,    1.0,  0,  1,  1,
  113.     kaiser_init, kaiser_print, (char *)&kd},
  114.    {0}
  115. };
  116.  
  117. /*-------------------- general filter routines --------------------*/
  118.  
  119. static filt_init()
  120. {
  121.     /* count the filters to set nfilt */
  122.     for (nfilt=0; nfilt<NFILTMAX && filt[nfilt].name; nfilt++);
  123.     /* better way?? */
  124.     mitchell_init(1./3., 1./3., (char *)&md);
  125.     kaiser_init(6.5, 0., (char *)&kd);
  126. }
  127.  
  128. /*
  129.  * filt_find: return ptr to filter descriptor given filter name
  130.  */
  131.  
  132. Filt *filt_find(name)
  133. char *name;
  134. {
  135.     int i;
  136.  
  137.     if (nfilt==0) filt_init();
  138.     for (i=0; i<nfilt; i++)
  139.     if (str_eq(name, filt[i].name))
  140.         return &filt[i];
  141.     return 0;
  142. }
  143.  
  144. /*
  145.  * filt_insert: insert filter f in filter collection
  146.  */
  147.  
  148. void filt_insert(f)
  149. Filt *f;
  150. {
  151.     if (nfilt==0) filt_init();
  152.     if (filt_find(f->name)!=0) {
  153.     fprintf(stderr, "filt_insert: there's already a filter called %s\n",
  154.         f->name);
  155.     return;
  156.     }
  157.     if (nfilt>=NFILTMAX) {
  158.     fprintf(stderr, "filt_insert: too many filters: %d\n", nfilt+1);
  159.     return;
  160.     }
  161.     filt[nfilt++] = *f;
  162. }
  163.  
  164. /*
  165.  * filt_catalog: print a filter catalog to stdout
  166.  */
  167.  
  168. void filt_catalog()
  169. {
  170.     int i;
  171.     Filt *f;
  172.  
  173.     if (nfilt==0) filt_init();
  174.     for (i=0; i<nfilt; i++)
  175.     filt_print(&filt[i]);
  176. }
  177.  
  178. /*
  179.  * filt_print: print info about a filter to stdout
  180.  */
  181.  
  182. void filt_print(f)
  183. Filt *f;
  184. {
  185.     fprintf(stderr,"%-9s\t%4.2f%s",
  186.     f->name, f->supp, f->windowme ? " (windowed by default)" : "");
  187.     if (f->printproc) {
  188.     fprintf(stderr,"\n    ");
  189.     filt_print_client(f);
  190.     }
  191.     fprintf(stderr,"\n");
  192. }
  193.  
  194. /*-------------------- windowing a filter --------------------*/
  195.  
  196. /*
  197.  * filt_window: given an IIR filter f and the name of a window function,
  198.  * create a compound filter that is the product of the two:
  199.  * wf->func(x) = f->func(x) * w->func(x/s)
  200.  *
  201.  * note: allocates memory that is (probably) never freed
  202.  */
  203.  
  204. Filt *filt_window(f, windowname)
  205. Filt *f;
  206. char *windowname;
  207. {
  208.     Filt *w, *wf;
  209.     window_data *d;
  210.  
  211.     if (str_eq(windowname, "box")) return f;    /* window with box is NOP */
  212.     w = filt_find(windowname);
  213.     ALLOC(wf, Filt, 1);
  214.     *wf = *f;
  215.     ALLOC(wf->name, char, 50);
  216.     sprintf(wf->name, "%s*%s", f->name, w->name);
  217.     wf->func = window_func;
  218.     wf->initproc = 0;
  219.     if (f->printproc || w->printproc) wf->printproc = window_print;
  220.     else wf->printproc = 0;
  221.     ALLOC(d, window_data, 1);
  222.     d->filter = f;
  223.     d->window = w;
  224.     wf->clientdata = (char *)d;
  225.     return wf;
  226. }
  227.  
  228. static double window_func(x, d)
  229. double x;
  230. char *d;
  231. {
  232.     register window_data *w;
  233.  
  234.     w = (window_data *)d;
  235. #   ifdef DEBUG
  236.     fprintf(stderr, "%s*%s(%g) = %g*%g = %g\n",
  237.     w->filter->name, w->window->name, x);
  238.     filt_func(w->filter, x), filt_func(w->window, x/w->filter->supp),
  239.     filt_func(w->filter, x) * filt_func(w->window, x/w->filter->supp));
  240. #   endif
  241.     return filt_func(w->filter, x) * filt_func(w->window, x/w->filter->supp);
  242. }
  243.  
  244. static void window_print(d)
  245. char *d;
  246. {
  247.     window_data *w;
  248.  
  249.     w = (window_data *)d;
  250.     if (w->filter->printproc) filt_print_client(w->filter);
  251.     if (w->window->printproc) filt_print_client(w->window);
  252. }
  253.  
  254. /*--------------- unit-area filters for unit-spaced samples ---------------*/
  255.  
  256. /* all filters centered on 0 */
  257.  
  258. double filt_box(x, d)        /* box, pulse, Fourier window, */
  259. double x;            /* 1st order (constant) b-spline */
  260. char *d;
  261. {
  262.     if (x<-.5) return 0.;
  263.     if (x<.5) return 1.;
  264.     return 0.;
  265. }
  266.  
  267. double filt_triangle(x, d)    /* triangle, Bartlett window, */
  268. double x;            /* 2nd order (linear) b-spline */
  269. char *d;
  270. {
  271.     if (x<-1.) return 0.;
  272.     if (x<0.) return 1.+x;
  273.     if (x<1.) return 1.-x;
  274.     return 0.;
  275. }
  276.  
  277. double filt_quadratic(x, d)    /* 3rd order (quadratic) b-spline */
  278. double x;
  279. char *d;
  280. {
  281.     double t;
  282.  
  283.     if (x<-1.5) return 0.;
  284.     if (x<-.5) {t = x+1.5; return .5*t*t;}
  285.     if (x<.5) return .75-x*x;
  286.     if (x<1.5) {t = x-1.5; return .5*t*t;}
  287.     return 0.;
  288. }
  289.  
  290. double filt_cubic(x, d)        /* 4th order (cubic) b-spline */
  291. double x;
  292. char *d;
  293. {
  294.     double t;
  295.  
  296.     if (x<-2.) return 0.;
  297.     if (x<-1.) {t = 2.+x; return t*t*t/6.;}
  298.     if (x<0.) return (4.+x*x*(-6.+x*-3.))/6.;
  299.     if (x<1.) return (4.+x*x*(-6.+x*3.))/6.;
  300.     if (x<2.) {t = 2.-x; return t*t*t/6.;}
  301.     return 0.;
  302. }
  303.  
  304. double filt_catrom(x, d)    /* Catmull-Rom spline, Overhauser spline */
  305. double x;
  306. char *d;
  307. {
  308.     if (x<-2.) return 0.;
  309.     if (x<-1.) return .5*(4.+x*(8.+x*(5.+x)));
  310.     if (x<0.) return .5*(2.+x*x*(-5.+x*-3.));
  311.     if (x<1.) return .5*(2.+x*x*(-5.+x*3.));
  312.     if (x<2.) return .5*(4.+x*(-8.+x*(5.-x)));
  313.     return 0.;
  314. }
  315.  
  316. double filt_gaussian(x, d)    /* Gaussian (infinite) */
  317. double x;
  318. char *d;
  319. {
  320.     return exp(-2.*x*x)*sqrt(2./PI);
  321. }
  322.  
  323. double filt_sinc(x, d)        /* Sinc, perfect lowpass filter (infinite) */
  324. double x;
  325. char *d;
  326. {
  327.     return x==0. ? 1. : sin(PI*x)/(PI*x);
  328. }
  329.  
  330. double filt_bessel(x, d)    /* Bessel (for circularly symm. 2-d filt, inf)*/
  331. double x;
  332. char *d;
  333. {
  334.     /*
  335.      * See Pratt "Digital Image Processing" p. 97 for Bessel functions
  336.      * zeros are at approx x=1.2197, 2.2331, 3.2383, 4.2411, 5.2428, 6.2439,
  337.      * 7.2448, 8.2454
  338.      */
  339.     return x==0. ? PI/4. : j1(PI*x)/(2.*x);
  340. }
  341.  
  342. /*-------------------- parameterized filters --------------------*/
  343.  
  344. double filt_mitchell(x, d)    /* Mitchell & Netravali's two-param cubic */
  345. double x;
  346. char *d;
  347. {
  348.     register mitchell_data *m;
  349.  
  350.     /*
  351.      * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
  352.      * SIGGRAPH 88
  353.      */
  354.     m = (mitchell_data *)d;
  355.     if (x<-2.) return 0.;
  356.     if (x<-1.) return m->q0-x*(m->q1-x*(m->q2-x*m->q3));
  357.     if (x<0.) return m->p0+x*x*(m->p2-x*m->p3);
  358.     if (x<1.) return m->p0+x*x*(m->p2+x*m->p3);
  359.     if (x<2.) return m->q0+x*(m->q1+x*(m->q2+x*m->q3));
  360.     return 0.;
  361. }
  362.  
  363. static void mitchell_init(b, c, d)
  364. double b, c;
  365. char *d;
  366. {
  367.     mitchell_data *m;
  368.  
  369.     m = (mitchell_data *)d;
  370.     m->p0 = (  6. -  2.*b        ) / 6.;
  371.     m->p2 = (-18. + 12.*b +  6.*c) / 6.;
  372.     m->p3 = ( 12. -  9.*b -  6.*c) / 6.;
  373.     m->q0 = (         8.*b + 24.*c) / 6.;
  374.     m->q1 = (      - 12.*b - 48.*c) / 6.;
  375.     m->q2 = (         6.*b + 30.*c) / 6.;
  376.     m->q3 = (     -     b -  6.*c) / 6.;
  377. }
  378.  
  379. static void mitchell_print(d)
  380. char *d;
  381. {
  382.     mitchell_data *m;
  383.  
  384.     m = (mitchell_data *)d;
  385.     fprintf(stderr,"mitchell: p0=%g p2=%g p3=%g q0=%g q1=%g q2=%g q3=%g\n",
  386.     m->p0, m->p2, m->p3, m->q0, m->q1, m->q2, m->q3);
  387. }
  388.  
  389. /*-------------------- window functions --------------------*/
  390.  
  391. double filt_hanning(x, d)    /* Hanning window */
  392. double x;
  393. char *d;
  394. {
  395.     return .5+.5*cos(PI*x);
  396. }
  397.  
  398. double filt_hamming(x, d)    /* Hamming window */
  399. double x;
  400. char *d;
  401. {
  402.     return .54+.46*cos(PI*x);
  403. }
  404.  
  405. double filt_blackman(x, d)    /* Blackman window */
  406. double x;
  407. char *d;
  408. {
  409.     return .42+.50*cos(PI*x)+.08*cos(2.*PI*x);
  410. }
  411.  
  412. /*-------------------- parameterized windows --------------------*/
  413.  
  414. double filt_kaiser(x, d)    /* parameterized Kaiser window */
  415. double x;
  416. char *d;
  417. {
  418.     /* from Oppenheim & Schafer, Hamming */
  419.     kaiser_data *k;
  420.  
  421.     k = (kaiser_data *)d;
  422.     return bessel_i0(k->a*sqrt(1.-x*x))*k->i0a;
  423. }
  424.  
  425. static void kaiser_init(a, b, d)
  426. double a, b;
  427. char *d;
  428. {
  429.     kaiser_data *k;
  430.  
  431.     k = (kaiser_data *)d;
  432.     k->a = a;
  433.     k->i0a = 1./bessel_i0(a);
  434. }
  435.  
  436. static void kaiser_print(d)
  437. char *d;
  438. {
  439.     kaiser_data *k;
  440.  
  441.     k = (kaiser_data *)d;
  442.     fprintf(stderr,"kaiser: a=%g i0a=%g\n", k->a, k->i0a);
  443. }
  444.  
  445. double bessel_i0(x)
  446. double x;
  447. {
  448.     /*
  449.      * modified zeroth order Bessel function of the first kind.
  450.      * Are there better ways to compute this than the power series?
  451.      */
  452.     register int i;
  453.     double sum, y, t;
  454.  
  455.     sum = 1.;
  456.     y = x*x/4.;
  457.     t = y;
  458.     for (i=2; t>EPSILON; i++) {
  459.     sum += t;
  460.     t *= (double)y/(i*i);
  461.     }
  462.     return sum;
  463. }
  464.  
  465. /*--------------- filters for non-unit spaced samples ---------------*/
  466.  
  467. double filt_normal(x, d)    /* normal distribution (infinite) */
  468. double x;
  469. char *d;
  470. {
  471.     /*
  472.      * normal distribution: has unit area, but it's not for unit spaced samples
  473.      * Normal(x) = Gaussian(x/2)/2
  474.      */
  475.     return exp(-x*x/2.)/sqrt(2.*PI);
  476. }
  477.